% 参数设置
alpha = 0.0001; % 热扩散率（常数）
T0 = 25; % 初始温度
c = 3.587; % 电极参数
P = 30; % 电极功率（W），假设恒定
Rmax = 20; % 网格最大半径
Zmax = 20; % 网格最大深度
dr = 0.1; % 网格步长（r方向）
dz = 0.1; % 网格步长（z方向）
dt = 0.1; % 时间步长
Tend = 8; % 模拟总时间

% 空间离散化
r = 0:dr:Rmax;
z = 0:dz:Zmax;
Nr = length(r);
Nz = length(z);

% 温度矩阵初始化
T = T0 * ones(Nr, Nz); % 初始温度设置为 T0

% 时间步数
nt = Tend / dt;

% 热源函数 Q(r, z, t)，假设电极功率在一定时间内持续不变
Q = @(r, z, t) (P / (c^3)); % 电极加热区域的功率密度（简化模型）

% 时间循环
for n = 1:nt
    T_new = T; % 创建新温度矩阵
    
    for i = 2:Nr-1
        for j = 2:Nz-1
            % 计算热源 Q(r, z, t) (假设恒定功率，持续加热)
            Q_val = Q(r(i), z(j), n*dt);
            
            % 使用有限差分法离散化热传导方程
            T_new(i,j) = T(i,j) + alpha * dt * ...
                ((T(i+1,j) - 2*T(i,j) + T(i-1,j)) / dr^2 + ...
                (T(i,j+1) - 2*T(i,j) + T(i,j-1)) / dz^2) + Q_val * dt;
        end
    end
    
    % 更新温度矩阵
    T = T_new;
    
    % 可视化温度分布（可选择性保存图像）
    if mod(n, 100) == 0 % 每10步保存一次图像
        figure;
        contourf(r, z, T', 20);
        colorbar;
        title(['Temperature Distribution at t = ', num2str(n*dt)]);
        xlabel('r');
        ylabel('z');
    end
end
